#include <TROOT.h>

void Z4F_CHIRALVAL(char *momentum_order, int beta100, int l, int t, int nclust, char *dir, int imom, char  *ensemble, char *pmprint, int useMPS, int subGP2)
{
	printf("=================  imom=%d  ===============\n",imom);
        L=l;
        T=t;
        Nclust=nclust;

        printf("Chiral fit for Dij");
        if (useMPS) printf(" using MPS^2");
        else printf(" using mpolar");
	printf("\n");

        double ZAPREV=0.757;
        int Nmass;



        // mass input files

        char mu_dat[1024];
        sprintf(mu_dat,"../INPUT/mu_%s%s_%d.dat",ensemble,pmprint,beta100);
        FILE *mudat=fopen(mu_dat,"r");
        fscanf(mudat,"%i",&Nmass);
        double *mu=(double*)malloc(sizeof(double)*(Nmass));
        for (int im=0; im<Nmass; im++)
        fscanf(mudat,"%lf",&mu[im]);


       char mPCAC_dat[1024];
       sprintf(mPCAC_dat,"../INPUT/mPCAC_%s%s_%d.dat", ensemble,pmprint,beta100);
       double mPCAC;
       FILE *mPCACdat=fopen(mPCAC_dat,"r");
       fscanf(mPCACdat,"%lf",&mPCAC);
       fscanf(mPCACdat,"%lf",&ZAPREV);


       if (useMPS)
       {
        int aux_int;
        double aux_double;
        char MPS_dat[1024];
        sprintf(MPS_dat,"../INPUT/MPS_%s%s_%d.dat",ensemble,pmprint,beta100);
        FILE *MPSdat=fopen(MPS_dat,"r");

        double **MPS=(double**)malloc(sizeof(double*)*(Nmass));
        for (int im1=0; im1<Nmass; im1++) MPS[im1]=(double*)malloc(sizeof(double)*(Nmass));

        for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++){
        fscanf(MPSdat,"%lf ",&MPS[im1][im2]);
        }
       }

      /////////////////////////
       //// Momenta set
      //////////////////////////////

    char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
    read_momenta_data(momentafile);
    double *nn=(double*)malloc(sizeof(double)*(Nmom));
    compute_nn(momentum_order, nn);
    Momentum(momentum_order);
    // equivalent momenta
    Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
    Average_over_Equivalent_Momenta(SinP2,SinP2,SinP2_eq);
 
 //////////////////////////////////
 // IO files
  ///////////////////////////////

  char file[1024];
  char dir_qav[1024];

 sprintf(dir_qav,"%s/QAV/",dir);

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

        //////////////////////////////////
        // Read Ds
       // I work with reals from now on
        /////////////////////////////
     double ****DpPC=malloc(sizeof(double***)*(5));
     for (int iop1=0; iop1<5; iop1++) DpPC[iop1]=malloc(sizeof(double**)*(5));
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DpPC[iop1][iop2]=malloc(sizeof(double*)*(Nmass));
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) for (int imass=0; imass<Nmass; imass++) DpPC[iop1][iop2][imass]=malloc(sizeof(double)*(Nclust+1));



     double *aux_clust=malloc(sizeof(double)*(Nclust+1));

     for (int imass=0; imass<Nmass; imass++)
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
     {
                sprintf(file,"DpPC%d%d_m%d",iop1+1,iop2+1,imass);
		//printf("%s/p%03d/%s\n",dir_qav,imom,file);
                read_input_onemomenta(dir_qav,file,aux_clust,imom);
                for (int iclust=0; iclust<Nclust+1; iclust++) DpPC[iop1][iop2][imass][iclust]=aux_clust[iclust];
     }


      ///////////////////////////////////////
      // CHIRAL FIT FOR EACH MOMENTUM 
     ///////////////////////////////////////

      double *ZAUX=malloc(sizeof(double)*(Nmass));
      double **ZAUX_CLUST=malloc(sizeof(double*)*(Nmass));
      for (int imass=0; imass<Nmass; imass++) ZAUX_CLUST[imass]=malloc(sizeof(double)*(Nclust+1));

      double *MAUX=malloc(sizeof(double)*(Nmass));

 


      double ***DCHI=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  DCHI[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DCHI[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));      

      double ***DLIN=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  DLIN[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DLIN[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

      double ***GP1=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  GP1[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) GP1[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

      double ***GP2=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  GP2[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) GP2[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

      double ***DCHI2=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  DCHI2[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DCHI2[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

	

     // chiral fit of D 
    // Substract simple pole for the D_1i and double for D_ij i =! 1
    // see DOCS/Frezzotti


    for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++){
	printf("D%d%d\n",iop1+1,iop2+1);
    for (int iclust=0; iclust<Nclust+1; iclust++)
     {
     for (int imass=0; imass<Nmass; imass++)
     {
		if (useMPS) MAUX[imass]=pow(MPS[imass][imass],2);
                else MAUX[imass]=sqrt(pow(mu[imass],2)+pow(ZAPREV*mPCAC,2));
                ZAUX[imass]=DpPC[iop1][iop2][imass][iclust];
                for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=DpPC[iop1][iop2][imass][iclust_loop];
     }
     if (iop1==0 || iop2==0 || ( iop1==1 && (iop2==3 || iop2==4)) || ( iop1==2 && (iop2==3 || iop2==4)) || ( iop1==3 && (iop2==1 || iop2==2)) || ( iop1==4 && (iop2==1 || iop2==2 ) )    )
     {
		if (iclust==Nclust) printf(" 1GP\n");
     GP2[iop1][iop2][iclust]=0.;
     lin_chiralfit_with1GP(ZAUX,ZAUX_CLUST,MAUX,Nmass,Nclust,DCHI[iop1][iop2][iclust],DLIN[iop1][iop2][iclust],GP1[iop1][iop2][iclust], DCHI2[iop1][iop2][iclust]); 
     sprintf(file,"DpPC%d%d",iop1+1,iop2+1); write_output_chiral_extrapolation_1GP(dir_qav,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass,DCHI[iop1][iop2],DLIN[iop1][iop2],GP1[iop1][iop2],DCHI2[iop1][iop2]);
      }

      else 
      {

	if (subGP2)
	{
		if (iclust==Nclust) printf(" 2GP\n");
	     lin_chiralfit_with2GP(ZAUX,ZAUX_CLUST,MAUX,Nmass,Nclust,DCHI[iop1][iop2][iclust],DLIN[iop1][iop2][iclust],GP1[iop1][iop2][iclust], GP2[iop1][iop2][iclust], DCHI2[iop1][iop2][iclust]);
	     sprintf(file,"DpPC%d%d",iop1+1,iop2+1); write_output_chiral_extrapolation_2GP(dir_qav,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass,DCHI[iop1][iop2],DLIN[iop1][iop2],GP1[iop1][iop2],GP2[iop1][iop2],DCHI2[iop1][iop2]);
	}
	else
	{
		if (iclust==Nclust) printf(" 1GP\n");
     	GP2[iop1][iop2][iclust]=0.;
     lin_chiralfit_with1GP(ZAUX,ZAUX_CLUST,MAUX,Nmass,Nclust,DCHI[iop1][iop2][iclust],DLIN[iop1][iop2][iclust],GP1[iop1][iop2][iclust], DCHI2[iop1][iop2][iclust]);
     sprintf(file,"DpPC%d%d",iop1+1,iop2+1); write_output_chiral_extrapolation_1GP(dir_qav,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass,DCHI[iop1][iop2],DLIN[iop1][iop2],GP1[iop1][iop2],DCHI2[iop1][iop2]);
	}
      }
     }
    }

    // write output
    printf("Writing output in %s/p%03d/(dat)(jk)\n",dir_qav,imom);

  //D_ij
   for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
   {
       	sprintf(file,"DpPC%d%d",iop1+1,iop2+1);
	write_output_onemomenta_real(dir_qav,file,DCHI[iop1][iop2],SinP2_eq,imom);
   }


    delete DpPC; 
    delete aux_clust;
    delete ZAUX; delete ZAUX_CLUST; delete MAUX;
    delete DCHI; delete DLIN; delete DCHI2; delete GP1; delete GP2;

}
